A First-Principles Method for Open Electronic Systems 
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We prove that the electron density function of a real physical system can be uniquely determined 
by its values on any finite subsystem. This establishes the existence of a rigorous density-functional 
theory for any open electronic system. By introducing a new density functional for dissipative 
interactions between the reduced system and its environment, we subsequently develop a time- 
dependent density-functional theory which depends in principle only on the electron density of the 
reduced system. In the steady-state limit, the conventional first-principles nonequilibrium Green's 
function formulation for the current is recovered. A practical scheme is proposed for the new density 
functional: the wide-band limit approximation, which is applied to simulate the transient current 
through a model molecular device. 



I. INTRODUCTION 

Density-functional theory (DFT) has been widely used 
as a research tool in condensed matter physics, chemistry, 
materials science, and nanoscience. The Hohenberg- 
Kohn theorem lays the foundation of DFT. The Kohn- 
Sham formalism 2] provides a practical solution to cal- 
culate the ground state properties of electronic systems. 
Runge and Gross extended DFT further to calculate the 
time-dependent properties and hence the excited state 
properties of any electronic systems • The accuracy of 
DFT or time-dependent DFT (TDDFT) is determined 
by the exchange-correlation (XC) functional. If the ex- 
act XC functional were known, the Kohn-Sham formal- 
ism would have provided the exact ground state prop- 
erties, and the Runge-Gross extension, TDDFT, would 
have yielded the exact time-dependent and excited states 
properties. Despite their wide range of applications, DFT 
and TDDFT have been mostly limited to isolated sys- 
tems. 

Many systems of current research interest are open 
systems. A molecular electronic device is one such sys- 
tem. DFT-based simulations have been carried out on 

such devices 0, 0, 0, S IM, 0, E3, EE These sim- 

ulations focus on steady-state currents under bias volt- 
ages. Two types of approaches have been adopted. One is 
the Lippmann-Schwinger formalism by Lang and cowork- 
ers Q. The other is the first-principles nonequilibrium 
Green's function (NEGF) technique @ |1 El El El • In 
both approaches the Kohn-Sham Fock operator is taken 
as the effective single-electron model Hamiltonian, and 
the transmission coefficients are calculated within the 
noninteracting electron model. The investigated systems 
are not in their ground states, and applying ground state 
DFT formalism for such systems is only an approxima- 
tion |13| . DFT formalisms adapted for current-carrying 
systems have also been proposed recently, such as Kosov's 
Kohn-Sham equations with direct current |l4j and Burke 
et a/.'s Kohn-Sham master equation including dissipation 



to phonons E3- However, practical implementation of 
these formalisms requires the electron density function 
of the entire system. In this paper, we present a rigor- 
ous DFT formalism for open electronic systems, and use 
it to simulate the steady and transient currents through 
molecular electronic devices. The first-principles formal- 
ism depends only on the electron density function of the 
reduced system. 

This paper is organized as follows. In Sec. |H]we pro- 
pose a TDDFT formalism for open electronic systems 
based on the equation of motion (EOM) for reduced 
single-electron density matrix. In Sec. 1 1 i 1 1 we prove the 
theorem that the electron density function of any finite 
subsystem can determine uniquely all properties of a con- 
nected real physical system. By utilizing this theorem we 
introduce in Sec. II VI a dissipation functional for the elec- 
tron density of the subsystem, and thus establish a rigor- 
ous and efficient first-principles formalism for steady and 
transient dynamics of open electronic systems. An wide- 
band limit (WBL) approximation scheme for the dissipa- 
tion functional is proposed for practical implementations 
in Sec. [V] To demonstrate the applicability of our first- 
principles formalism, a TDDFT calculation is carried out 
to simulate the transient current through a model molec- 
ular device. The detailed procedures and results are de- 
scribed in Sec. IVII Discussion and summary are given in 
Sec.rvTil 
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II. FIRST-PRINCIPLES FORMALISM 

A. Reduced single-electron density matrix and 
TDDFT formalism for reduced system 

Fig- LH depicts an open electronic system. Region D is 
the reduced system of our interests, and the electrodes 
L and R are the environment. Altogether D, L and R 
form the entire system. Taking Fig. ^ as an example, 
we develop a practical DFT formalism for the open sys- 
tems. Within the TDDFT formalism, a closed EOM has 
been derived for the reduced single-electron density ma- 
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trix a(t) of the entire system [l6j : 

i&(t) = [h(t),a{t)], 



(1) 



where h(t) is the Kohn-Sham Fock matrix of the entire 
system, and the square bracket on the right-hand side 
(RHS) denotes a commutator. The matrix element of a 
is defined (t) = (aUt) cii(t)}, where a,i(t) and a\(t) 

are the annihilation and creation operators for atomic 
orbitals i and j at time t, respectively. Fourier trans- 
formed into frequency domain while considering linear 
response only, Eq. Q leads to the conventional Casida's 
equation [13. Expanded in the atomic orbital basis set, 
the matrix representation of a can be partitioned as 



&l old olr 

&DL &D &DR 

&RL ord or 



(2) 



where <jl, <jr and <td represent the diagonal blocks cor- 
responding to the left lead L, the right lead R and the de- 
vice region D, respectively; old is the off-diagonal block 
between L and D; and ord, olr, <Jdl, &dr. and grl 
are similarly defined. The Kohn-Sham Fock matrix h 
can be partitioned in the same way with a replaced by h 
in Eq. (|2J). Thus, the EOM for ctd can be written as 



io D 



[hD,o D ] 
[hD,cr D ] 



E 

a=L,R 



{h Da G aD — <J Da h aD ) 
E Oo. ( 3 ) 



a=L,R 



where Ql (Qr) is the dissipative term due to L (i?). 
With the reduced system D and the leads L/R spanned 
respectively by atomic orbitals {1} and single-electron 
states {k a }, Eq. J2J is equivalent to: 



leD 

hot £a 



{hnlOlm — Onlhlm) 



t ^ ^ Qa,n7n> (4) 
a=L,R 



tk n Ok n 



Onk„ 



0> 



(5) 



where m and n correspond to the atomic orbitals in re- 
gion D; k a corresponds to an electronic state in the elec- 
trode a (a = L or R). h nka is the coupling matrix el- 
ement between the atomic orbital n and the electronic 
state k a . The current through the interfaces Sl or Sr 
(see Fig. ^) can be evaluated as follows, 



J a (t) 



, 9 , x 
dr-QiP(r,t) 



E | **<.*.(*) 

k a ea 



EE 

leD k a £c 

= -^Q a ,u = -tr[Q a (t)] 



- > > [h ka i <Ji ka 

l&D k a £a 



Ok a l h 



lk, x 



(6) 



leD 

i.e., the trace of Q a . 



B. Solution for steady-state current 

Based on the Keldysh formalism and the analytic 
continuation rules of Langreth 0], Q a ,nm{t) can be cal- 
culated by the NEGF formulation as described in Refer- 
ence [2(J (see Appendix |A"|) 



-E/ H 

leD J-oo 

G r nl (t,T)X< lm (T 7 t) - X< nl (t,T)G? m (T,t) 



(7) 



where G r , G a and G < are the retarded, advanced and 
lesser Green's function for the reduced system D, respec- 
tively, and S^,, and are the retarded, advanced 
and lesser self-energies due to the lead a (L or R), re- 
spectively. Combining Eqs. © and JJJ), we obtain 



J a {t) = 23? M drti 

I J — OO 

GJ,(*,r)£<(r,t) 



G<(t,r)SS(r,t) + 



(8) 



The same expression was derived by Stefanucci and Alm- 
bladh within the framework of TDDFT 

It is important to point out that Eqs. Q — GJ follow 
the partition-free scheme proposed by Cini |22j, while 
Eq. was derived if one follows the partitioned scheme 
developed by Caroli et al. |23j . In the above derivation we 
assume that the equivalence of the two schemes, which 
is satisfied if the two self-energies behave asymptotically 
as follows 211, 



Jim K(t,t') = Jim K(t',t) = 0. 



(9) 



As t, t -> +oo, T k n °; n (t,T) = h nka (t) h kam (r) becomes 
asymptotically time-independent. The Green's functions 
for the reduced system D rely simply on the difference of 
the two time- variables [2lJ , and can thus be expressed as 



g< 



(*,T) 



E 

p,q£D 



dt 1 



OO 



dt% G r np {t, ti) 



^< q {txM)G a qm {t 2 ,r) 

* E E E /r 

p.qeD a—L,R l a £a 
oo 

dtie-^G^-ft-ti) 



x 



OO 



dt 2 e ie ' i2 G« (t 2 -r) 



f E EE// 

p,q£D a=L,R I„ea 

x Gn„(ef ) r'° G a (e") 



-ie["(t-T) 



Jnm ' 



G5»(e) = [cJ-fc -E^c)] 

= E E^^-^iiA)- 1 



(10) 

(11) 
(12) 



a=L,R lea 



3 



where / is the identity matrix. The steady-state cur- 
rent can thus be explicitly expressed by combining 
Eqs. (ElJ-jnj), 



J L {co) 



-Jr(oo) 

- QL,nn(oo) 



n£D 



l k£L l£R 



x tr 



leR k£L 

r { r^i? r* a f r^x 



T(e) 



G r D (ef)r l «G a D (ef)r K 

[/ L (e)- f R (e)]T(e)de, (13) 
2 7 r ?7i ^tr[G^(e)r«(e)G^(e)r i (e)l.(14) 



Here T(e) is the transmission coefficient, f a (e) is the 
Fermi distribution function, and i] a (e) = J2 kea S(e — e£) 
is the density of states (DOS) for the lead a (L or 
R). Eq. (|T3*)) is exactly the Landauer formula HH 
in the DFT-NEGF formalism 0, E3- The difficulty in 
solving Eq. is to calculate Q a ,nm- Employing the 
Keldysh NEGF formalism, the evaluation of Q a ,nm in- 
volves the calculation of two-time Green's functions and 
self-energies as those appearing in Eq. Q, which makes 
the simulation of any real molecular device computation- 
ally impractical. An alternative approach must be devel- 
oped. 



III. HOLOGRAPHIC ELECTRON DENSITY 
THEOREM FOR TIME-DEPENDENT SYSTEMS 

As early as in 1981, Riess and Munch discovered 
the holographic electron density theorem which states 
that any nonzero volume piece of the ground state elec- 
tron density determines the electron density of a molec- 
ular system. This is based on that the electron den- 
sity functions of atomic and molecular eigenfunctions are 
real analytic away from nuclei. In 1999 Mezey extended 
the holographic electron density theorem |2g. And in 
2004 Fournais et al. proved again the real analyticity of 
the electron density functions of any atomic or molecular 
eigenstates [2!j. Therefore, for a time-independent real 
physical system made of atoms and molecules, its elec- 
tron density function is real analytic (except at nuclei) 
when the system is in its ground state, any of its excited 
eigenstates, or any state which is a linear combination 
of finite number of its eigenstates; and the ground state 
electron density on any finite subsystem determines com- 
pletely the electronic properties of the entire system. 



Molecular 
Device 



Switch 
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FIG. 1: Schematic representation of the experimental setup 
for quantum transport through a molecular device. 



As for time-dependent systems, the issue is less clear. 
Although it seems intuitive that the electron density 
function of any time-dependent real physical system is 
real analytic (except for isolated points in space-time), 
it turns out quite difficult to prove the analyticity rig- 
orously. Fortunately we are able to establish a one-to- 
one correspondence between the electron density function 
of any finite subsystem and the external potential field 
which is real analytic in both i-space and r-space, and 
thus circumvent the difficulty concerning the analyticity 
of time-dependent electron density function. For time- 
dependent real physical systems, we have the following 
theorem: 

Theorem: If the electron density function of a real 
physical system at to, p(r,to), is real analytic in r- 
space, the corresponding wave function is $(io)j aud the 
system is subjected to a real analytic (in both t-space 
and r-space) external potential field v(r,t), the time- 
dependent electron density function on any finite sub- 
space D, p£>(r,t), has a one-to-one correspondence with 
v(r, t) and determines uniquely all electronic properties 
of the entire time-dependent system. 

Proof: Let v(r, t) and v'(r,t) be two real analytic po- 
tentials in both i-space and r-space which differ by more 
than a constant at any time t ^ to, and their correspond- 
ing electron density functions are p(v,t) and p'(r,t), re- 
spectively. Therefore, there exists a minimal nonnegative 
integer k such that the fc-th order derivative differentiates 
these two potentials at to- 



±L[v(r,t)-v'(r,t)} 



^ const. 



(15) 



t=to 



Following exactly the Eqs. (3)-(6) of Ref. |{g, we have 



Qk+2 

dt k + 2 



[p(r,t)-p'(r,t)] 



-V-it(r), (16) 



t=t 



where 



u(r) = p(r,t )V< 



[v(r,t)-v'(r,t)] 



.(17) 



Due to the analyticity of p(r,t ), v(r, t) and v'(r,t), 
V-it(r) is also real analytic in r-space. It has been proven 
in Ref. Q that it is impossible to have V ■ u(r) = on 
the entire r-space. Therefore it is also impossible that 
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V • u(r) — everywhere in D because of analytical con- 
tinuation of V • u(r). Note that po{r,t) = p(r,t) for 
re J). We have thus 



s 



Qk+2 
dt k + 2 



[p D (r,t) - p' D (r,t)} 



+ 



(18) 



for r G D. This confirms the existence of a one-to- 
one correspondence between v(r, t) and p_o(r, t). /0£>(r, t) 
thus determines uniquely all electronic properties of the 
entire system. This completes the proof of the Theorem. 

Note that if <E»(£ ) is the ground state, any excited 
eigenstate, or any state as a linear combination of fi- 
nite number of eigenstates of a time-independent Hamil- 
tonian, the prerequisite condition in Theorem that the 
electron density function p(r, to) be real analytic is auto- 
matically satisfied, as proven in Ref. 29]. As long as the 
electron density function at t = to, p(r,to), is real ana- 
lytic, it is guaranteed that pjj(r,t) of the subsystem D 
determines all physical properties of the entire system at 
any time t if the external potential v(r, t) is real analytic. 

According to the above Theorem, the electron density 
function of any subsystem determines all the electronic 
properties of the entire time-dependent physical system. 
This proves in principle the existence of a rigorous DFT- 
type formalism for open electronic systems. All one needs 
to know is the electron density of the reduced system. 



IV. DISSIPATIVE DENSITY FUNCTIONAL 

According to the holographic electron density theorem 
of time-dependent physical systems, all physical quan- 
tities are explicit or implicit functionals of the electron 
density in the reduced system D, p D (r,t). Q a of Eq. © 
is thus also a universal functional of pu(r,t). Therefore, 
Eq. J2J can be recast into a formally closed form, 



ia D 



h D [r,t;p D (r,t)],a D - i Q a [r , t; p D (r, t)\ . 



-L.R 



(19) 

Neglecting the second term on the RHS of Eq. (|19fl leads 
to the conventional TDDFT formulation in terms of re- 
duced single-electron density matrix 0] for the isolated 
reduced system. The second term describes the dissi- 
pative processes between D and L or R. Besides the 
XC functional, an additional universal density functional, 
the dissipation functional Q a [r, i;p,o(r, t)], is introduced 
to account for the dissipative interaction between the 
reduced system and its environment. Eq. I|19|) is the 
TDDFT EOM for open electronic systems. It would thus 
be much more efficient integrating Eq. (|19f) than solving 
Eqs. and J7J, if Q a [r,t; pr>(r,t)] or its approxima- 
tion is known. We therefore have a practical and poten- 
tially rigorous formalism for any open electronic systems. 
Burke et al. extended TDDFT to include electronic sys- 
tems interacting with phonon baths |15| , they proved the 
existence of a one-to-one correspondence between v(r, t) 
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FIG. 2: Schematic representation of the reduced system D 
within the WBL scheme for dissipation functional Q a . 



and p(r, t) under the condition that the dissipative inter- 
actions (denoted by a superoperator C in Ref. 0]) be- 
tween electrons and phonons are fixed. In our case since 
the electrons can move in and out the reduced system, 
the number of the electrons in the reduced system is not 
conserved. In addition, the dissipative interactions can 
be determined in principle by the electron density of the 
reduced system. We do not need to stipulate that the 
dissipative interactions with the environment are fixed 
as Burke et al.. And the only information we need is 
the electron density of the reduced system. In the frozen 
DFT approach [3(| an additional XC functional term was 
introduced to account for the XC interaction between the 
system and the environment. This additional term is in- 
cluded in h D [r,t; p D (r,t)] of Eq. (|T§|| . 

Given Q a [p D (r,t)\ how do we solve the EOM {T^Jl 
in practice? Again take the molecular device shown in 
Fig. ^ as an example. We focus on the reduced system 
D as depicted in Fig. |2 and integrate the EOM (fCT?)l di- 
rectly by satisfying the boundary conditions at Sl and 
Sr. AV L (t) and AV R (t) are the bias voltages applied 
on L and R, respectively, and serve as the boundary 
conditions at Sl and Sr, respectively. At t — > — oo, 
AV L = AV R = 0, and near t = AV L (t) and AV R (t) 
are turned on adiabatically. We need thus integrate 
Eq. (|19fl together with a Poisson equation for Coulomb 
potential inside the device region D. And the Poisson 
equation is subjected to the boundary condition deter- 
mined by the potentials at Sl and Sr. 



V. WIDE-BAND LIMIT APPROXIMATION 
FOR DISSIPATION FUNCTIONAL Q a AND ITS 
TEST ON A MODEL SYSTEM 

An explicit form for the dissipation functional Q a is 
required for practical implementation of Eq. (|19|) . Ad- 
mittedly Q Q [r, t; pd{y, <)] is an extremely complex func- 
tional and difficult to evaluate. As various approximated 
expressions have been adopted for the DFT XC func- 
tional in practical implementations, the same strategy 
can be applied to the dissipation functional Q a . 

One such scheme is the wide-band limit (WBL) ap- 
proximation p(il ] which involves the following assump- 



tions for the leads: (i) their band-widths are assumed to 
be infinitely large, (ii) their line-widths, A£(i, t), defined 
by the DOS at Sl or Sr times the coupling strength be- 
tween D and L or R, i.e., A£(*,t) = tt r] a (e%) T k " (t, r), 
are treated as energy independent, i.e., A^(t,r) « 
A Q (i,r) « A Q , and (iii) the level shifts of L or i? are 
taken as a constant for all energy levels, i.e., Ae£(i) « 
Ae"(t) = -AV a (t), where AF Q (t) are bias voltages ap- 
plied on L or R at time f. The detailed derivations for 
the WBL scheme can be found in Appendix [B] an d the 
explicit form for Q^ BL is given here, 



JlO) 



Jn(t) 



Q^ Bi ®=lf°(t) + {A<> D (i)}, 
Here K a (t) is fully expanded as follows, 

dee iet 



(20) 



+ / [/-£/"(*)« 



where 



- ft D (t) +iA + Ae«(i) 



A«+ff.C(21) 



(22) 



From Eqs. I|20|) - (|22[l it is clear that the dissipation func- 
tional Q a within WBL scheme depends explicitly on 
Ae a (t), a D (t), h D {t) and A a . Note that Ae a {t) = 
-AV a (t) and AV a (t) is a functional of p D {r,t), i.e., 
AV a (t) = AV a [ PD (r,t),t]; h D (t) = h D [a D (t),t]; since 
77a (e) is the DOS at S a (a — L or R), and T Q is the cou- 
pling strength between the surface states at S a and the 
bulk states of D, A" is thus a functional of pu(r,t), i.e., 
A a = A a [pd (r, t),t]. We hence conclude that in practice 
qWBL - 1S a f unc tional of /Ou(r, i), i.e., 



iWBLi 



(t) = cT D [p D ],h D [a D [ PD },t], 



A a [p D ,t], AV a [ PD ,t], t 



(23) 



The WBL dissipation functional Q^f BL is then tested 
by calculations on a model system which has previously 
been investigated by Maciejko, Wang and Guo j33|. In 
this model system the device region D consists of a sin- 
gle site spanned by only one atomic orbital (see Fig.|3J). 
Exact transient current driven by a step voltage pulse 
has been obtained from NEGF simulations [3l| , and the 
authors concluded that the WBL approximation yields 
reasonable results provided that the band-widths of the 
leads are five times or larger than the coupling strength 
between D and L or R. The computational details are 
as follows. The entire system (L + R + D) is initially 
in its ground state with the chemical potential p°. Ex- 
ternal bias voltages are switched on from the time t = 0, 
which results in transient current flows through the leads 
L and R. 5h D {t) = h D {t) - h D (0), Ae L (t) and Ae R (t) 
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FIG. 3: Model system for the test of the WBL dissipation 
functional where a single site spans the device region D. Tran- 
sient currents through leads L and R, Jl(£) and Jn.it), are 
simulated. The inset shows the time-dependent level shift of 
lead R. 
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FIG. 4: The calculated transient current through Sr within 
the WBL scheme. We set p° = hn(0) = for the ground 
state; and Ae L (t) = and Ae R (t) = Ae R (1 - e~ t/a ) after 
switch-on. The above panels show different cases where (a) 
Ae R = 2 eV, A L = A R = 0.1 eV; (b) Ae R = 0.2 eV, A L = 
A R = 0.1 eV; (c) Ae R = 10 eV, A L = A R = 0.1 eV; and (d) 
Ae R = 2 eV, A L = A R = 0.04 eV, respectively. 



are the level shifts of D, L and R at time t, respectively. 
In our works we take Sh D (t) = \ [Ae L (t) + Ae R (t)] , 
Ae L (t) = 0, and Ae R (t) = Ae R (1 - e"*/ ), where a is 
a positive constant. The real analytic level shift Ae R (t) 
resembles perfectly a step pulse as a —> + . The calcu- 
lation results are demonstrated in Fig. 0] We choose ex- 
actly the same parameter set as that adopted for Fig. 2 in 
Ref. [3l] , and the resulting transient current, represented 
by Fig^IJa), excellently reproduces the WBL result in 
Ref. [31 1 , although the numerical procedures employed 
are distinctively different. The comparison confirms ev- 
idently the accuracy of our formalism. From Fig. 0Ja)- 
(c) it is observed that with the same line- widths A", a 
larger level shift Ae R results in a more fluctuating cur- 
rent, whereas by comparing (a) and (d) we see that un- 
der the same Ae R , the current decays more rapidly to 
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FIG. 5: A two-dimensional model molecular device is con- 
nected to left and right leads. 



the steady state value with larger A Q . 

Since the integration over energy in Eq. I|21[l can be 
performed readily by transforming the integrand into di- 
agonal representation, Q^ BL are evaluated efficiently, 
which makes the WBL scheme a practical routine for 
subsequent TDDFT calculations. 



VI. A TDDFT CALCULATION OF TRANSIENT 
CURRENT 



With the EOM (JTJJJl and the WBL scheme for the dis- 
sipation functional Q a , it is now straightforward to carry 
out first-principles calculations for transient dynamics of 
open electronic systems. A model molecular device de- 
picted in Fig. [3 is taken as the open system under in- 
vestigation. The device region D containing 24 carbon 
and 12 hydrogen atoms is spanned by the 6-31 Gaus- 
sian basis set, i.e., altogether 240 basis functions for the 
reduced system. The leads are quasi-one-dimensional 
graphene sheets with dangling bonds saturated by hy- 
drogen atoms, and the entire system is on a same plane. 
The ground state reduced single-electron density matrix 
for the reduced system, <td(0), is extracted from cr(0) of 
an extended system which consists of totally 134 atoms, 
covering not only the device region D but also portions of 
leads L and R. This provides the initial condition for the 
EOM O. The line-widths A L and A R within the WBL 
scheme are obtained from the surface Green's functions 
for isolated semi- infinite bulk leads L and R, g£(/i°) and 
g r R (n°) H3, respectively, and then optimized such that 
the RHS of the EOM vanishes correctly at t = 0. 
The molecular device is switched on by a step-like volt- 
age AV R (t) = -Ae R (t) = AV R (l-e-^ a ) applied on the 
right lead (see the inset of Fig-EJ) , while AV L {t) = 0, and 
the dynamic response of the reduced system is obtained 
by solving the EOM (|19fl in time domain within the adi- 
abatic local density approximation (ALDA) 01 f° r the 
XC functional and the WBL approximation for the dis- 
sipation functional. To save computational resources we 
linearize the XC component of the induced Kohn-Sham 
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FIG. 6: The solid (dotted) line represents the transient cur- 
rent through Sr (Sl) driven by a step-like voltage AV R (t) 
applied on the lead R with the amplitude (a) AV R = — 1 mV, 
and (b) AV R = -1 V. 



Fock matrix on D, 5hp C (i), as follows, 



XC 



E V ^'n Km(*) - <W(0)] , (24) 
Sv xc [p D (r,t)} 



mn£D 

dr^(r)^(r) 

D 

x#(r)&(r), 



dp D (r,t) 



(25) 



where v xc [p£>(r, t)] is the XC potential. The Coulomb 
component of Shz)(t) is constructed by solving the Pois- 
son equation for the device region D subjected to bound- 
ary conditions AV a (t) at every time t. The TDDFT 
calculations are carried out with a modified version of 
the TDDFT-LDM program developed by Yam, Yokojima 
and Chen UH. 

In Fig. HJa) and (b) we plot the transient currents 
through the interfaces Sl and Sr, Jh(t) and JRif), 
for cases where the turn-on voltage AV R = —1 mV 
and AV R = -1 V, respectively. The EOM O is in- 
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tegrated numerically by the fourth-order Runge-Kutta 
method [33] up to 25 fs with the time step 0.02 fs. 
Jh{t) and Ji?(t) depicted in Fig. UJa) increase rapidly 
in the first 5 fs and then approach gradually towards 
their steady state values. The steady currents through 
the leads L and R are —62.8 nA and 62.8 nA, respec- 
tively and thus cancel each other out exactly, as they 
should. With a much larger turn-on voltage AV R , Ji(t) 
and Jn(t) exhibit conspicuous overshooting during the 
first 2 fs, as shown in Fig. |5fb), and afterwards they de- 
cay slowly to their steady state values, i.e., —21.4 /iA 
and 21.4 ^A, respectively. From the both cases shown in 
Fig. diversified fluctuations are observed for the time- 
dependent currents. This is due to the various eigen- 
values possessed by the non-negative definite line-widths 
A™ with their magnitudes ranging from to 39 eV, corre- 
sponding to various dissipative channels between D and 
L or R. For much higher turn-on voltages the linearized 
form for 5hp(t) (Eq. (|24|l ) becomes inadequate, which 
makes such a TDDFT calculation computationally de- 
manding with our present coding. From Fig. EI the char- 
acteristic switch-on time for the model molecular device 
depicted in Fig. [S] is estimated as about 10 fs for applied 
bias voltages as large as 1 V. 



VII. DISCUSSION AND SUMMARY 

With an explicit form of the universal dissipation func- 
tional Q a , the time evolution of an open electron sys- 
tem in external fields is fully characterized by the EOM 
for the reduced single-electron density matrix of the re- 
duced system (see Eq. $Tty ). In practical calculations, 
we need thus focus only on the reduced system with ap- 
propriate boundary conditions. In conventional quantum 
dissipation theory (QDT) [34| the key quantity is the re- 
duced system density matrix. Whereas in Eq. (|19|l the 
basic variable is the reduced single-electron density ma- 
trix, which leads to the drastic reduction of the degrees of 
freedom in numerical simulation. Linear-scaling methods 
such as the localized-density- matrix method |l6ll35j | may 
thus be adopted to further speed up the solution process 
of Eq. (|19fl . Yokojima et al. developed a dynamic mean- 
field theory for dissipative interacting many-electron sys- 
tems [3(| |37} • An EOM for the reduced single-electron 
density matrix was derived to simulate the excitation and 
nonradiative relaxation of a molecule embedded in a ther- 
mal bath. This is in analogy to our case although our 
environment is actually a fermion bath instead of a bo- 
son bath. More importantly the number of electrons in 
the reduced system is conserved in Refs. [3^, |3tJ while in 
our case it is not. Therefore, Eq. Ijl9(l provides a rigor- 
ous and convenient formalism to investigate the dynamic 
properties of open systems. Recently Cui et al. pro- 
posed a TDDFT scheme for first-principles study of non- 
equilibrium quantum transport based on the complete 
second-order quantum dissipation theory (CS-QDT) [38| . 
their formulation is constructed in terms of an improved 



reduced density matrix approach at the self-consistent 
Born approximation (SCBA) level. 

It is worth mentioning that our first-principles method 
for open systems applies to the same phenomena, prop- 
erties or systems as those intended by Hohenberg and 
Kohn pj, Kohn and Sham Q, and Runge and Gross 0], 
i.e., where the exchange-correlation energy is a functional 
of electron density only Exc = Exc[p{ r )]- This is true 
when the interaction between the electric current and 
magnetic field is negligible. However, in the presence 
of a strong magnetic field, Exc = Exc[p( r ) Jp( r )] or 
Exc = Exc[p{ r ), B(r)], where j p (r) is the paramag- 
netic current density and B(r) is the magnetic field [39(. 
In such a case, our first-principles formalism needs to be 
generalized to include j p (r) or B(r). Of course, j p (r) or 
B(r) should be an analytical function in space. It is im- 
portant to note that our formalism applies in principle 
to Cini's scheme. Caroli's scheme is employed to derive 
an approximated expression for the dissipative functional 

Qa ■ 

To summarize, we have proved rigorously the existence 
of a first-principles method for time-dependent open elec- 
tronic systems, and developed a formally closed TDDFT 
formalism by introducing a new dissipation functional. 
This new functional Q a depends only on the electron 
density function of the reduced system. With an efficient 
WBL scheme for Q a , we have applied the first-principles 
formalism to carry out a TDDFT calculation for tran- 
sient current through a model molecular device. This 
work greatly extends the realm of density-functional the- 
ory. 
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APPENDIX A: DERIVATION OF EQ. © 

In Keldysh formalism 0]> the nonequilibrium single- 
electron Green's function Gk a , m {t,t') is defined by 

G kam (t,t') = -i(T c {a ka (t)al{t')}) , (Al) 

where Tq is the contour-ordering operator along the 
Keldysh contour [l^. Its lesser component, Gj5 m (t,t'), 
is defined by 

Gt m (t,t')=i{al(t')a ka (t)). (A2) 

Therefore <Jk a m{t) is precisely the lesser Green's func- 
tion of identical time variables, i.e., CFk a m(t) = 
-iG^ m (t, t')|t'=t- The formal NEGF theory has exactly 
the same structure as that of the time-ordered Green's 
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function at zero temperature. Thus, the Dyson equation 
for G kam (t,t') can be written as 

Gk a m(t, t') = / dr gk a {t,T)h ka i(T)Gi m {T,t'), 
leD Jc 

(A3) 

where Gj m (r, t') and g ka (t,r) are the contour-ordered 
Green's functions for the reduced system D and the iso- 
lated lead a (L or R), respectively. Applying the analytic £i fa £\ 
continuation rules of Langreth [lij , we have 

GLmW) = E / dTh kal (r) gt(t,r)GUr,t') 



leD 



+ 9L(t,T)Gf m (T,t') 



(A4) 



where Gf m (T, t') and Gf m {T,t') are the advanced and 
lesser Green's functions for the reduced system D, re- 
spectively, and g ka (t, r) and g k (t, r) are the retarded 
and lesser Green's functions for the isolated lead a (L or 
R) [2(j, respectively. Note that 

G< ka (t',t) = i(a{ a (t)a m (t')) 

= -[G< am (t,t')}\ (A5) 

Obviously a mka (t) = -iG< ka (t' ',t)|f=t. Combining 
Eqs. (|A4|) and (|A5|> . we obtain 

/>oo 



by employing the following equalities: 



(A6) 



G r ml (t',r) = [G? m (r,t')}* , 

GUt',r) = -[Gf m {T,t')]\ 

9t(r,t) = [g r ka (t,r)Y, 

9t(r,t) = -[gt(t,r)]\ 



(A7) 



By inserting Eqs. (|A^|l and HA6|> into Eq. (gj, Eq. J7J) is 

recovered straightforwardly where the self-energy terms 
are defined by 



E hkJt) 9kJ^ T ) ' h k a m{T), 

-a Get 

E h nk a (t)g ka (t,T)h ka i(T), 

E ^.WsKi.^wW- (A8) 



APPENDIX B: WIDE-BAND LIMIT 
APPROXIMATION FOR DISSIPATION 
FUNCTIONAL Q Q 



Within the WBL scheme, the retarded and advanced 
self-energies become local in time po| . 



h nk a {t)h kam (T)g1: a (T,t) 

fc Q 6ct 

^ h nka (t)h kam (T) 

k a Get 



+ 00 



de > A 



,(r,t) 



i^-r)A« m , 

[^a,mii(^ T )] 

- r)A« rr , 



(Bl) 
(B2) 



The third equality of Eq. (|B1|) involves the following ap- 
proximation for the line-widths within the WBL scheme, 



A£(i,r) = irr] a (e k )h nka (t)h kam (T) 

« A"(t,r) ~ A Q . (B3) 



Initially the entire system (D + L + R) is in its ground 
state with the chemical potential from the time t = 
it is switched on by external potentials AV a (t) applied 
on the leads L or R. Hence, for t, t > we have 



KkA^k^ii^g^ij.t) 

^ hnk a (t)h kam (T) 

= a 



k a ea 
X 



?l e »/* Ae^CQdt A a 



+ OC 



/a( e ) e fc (*- T )deL (B4) 



G r nm (t,r) = -i^t-T)J2Ki } (t)U^(r), (B5) 



k a £a 



where Ae a (t) = — AV a (t) are the time-dependent level 
shifts for the leads L and R, while for r < and t > 0, 
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the counterparts of (|B4(1 and (|B5(1 are as the following: 

Z< nm (r,t) = h nk a (t)h kam {r)g^(T,t) 



fc Q Set 



^ hnk a (t)hk a m(T) 



X 

2? 



pi/o Ae Q (f)dt 

1 'ran 



(B6) 



G r nm (t,r) = ^Utf® G{ m (0,T) 



(B7) 



where £?['(— r) is the retarded Green's function for the 
reduced system D before switch-on. The propagators for 
the reduced system U^'(t) are defined as 



expj ±ij h D (r)dT± Atj. (B8) 



where A = J2 a =L,n Aa - B V inserting Eqs. (jBT) - (|B7|) 
into Eq. (JJJ the explicit form of WBL approximation for 
the dissipation functional Q a is obtained as 

Q™ BL {t) = K a (t) + {A a ,a D (t)}, (B9) 

where the curly bracket on the RHS denotes an anticom- 
mutator, and K a {t) is a Hermitian matrix expressed by 



K a (t) = P a {t) + [P a (t)}'< , 



(BIO) 



where P a (t) involve an integration over the entire t- 
space, which is then decomposed into positive and nega- 
tive parts, denoted by Pi +S> (t) and P&(t), respectively. 

r+oo 

P°(t) = / dTG r D (t,T)X<(T,t) 



(Bll) 



Pt\t) and P<T ; (t) are evaluated via 



>(+) 



drG r D (t,T)E<(r,t) 



2i 



exp 



i Ae a (T)dT}U^-\t) 



e - h D (0) + iA 



A Q , (B12) 



and 



^0 

P«(t) = -- [ deWt\e,t) 



x / drW^ +) (e,T)A tt , (B13) 
Jo 



respectively, where 

W ± (e,t) = e ±i I° *>"[/lr 1 (T)-iA-Ae<*(<r)-e]_ (B14) 

However, the evaluations of Eqs. I)B13(I - I|B14|I are found 
extremely time-consuming since at every time t one needs 
to propagate W^(e,t) for every individual e inside the 
lead energy spectrum. It is thus necessary to seek for a 
simpler approximate form for Pi + \t) with satisfactory 
accuracy retained. Note that Eq. I|B13(1 can be reformu- 
lated as 

Pi +) (t) = -II del dr 



" J-oo JO 
x g-i/'l>.D(t)-*A-Ae '(t)-e]<iF^a > CQl^ 

For cases where a steady state can be ultimately reached, 
Ae a (t) and become asymptotically constant as 

time t -» +oo, i.e., Ae a (t) -> Ae Q (oo) and -> 
h D (oo). Therefore, the steady state P<i + ^(oo) can be 
approximated by substituting Ae Q (oo) and /id(co) for 
Ae"(i) and /loft) in Eq. i|ET5|) . respectively. 



oo 



x g— i[fcn(oo)— SA— Ae°(oo)- e](t— t) 

_2i r»° 



X 



_f y 'i[/ir> (oo) — iA— Ae" (oo) — e]t \ 

A". (B16) 



e - h D (oo) + iA + Ae a (oo) 

It is obvious from Eq. I|B15J| that 

Pi +) (0)=0. (B17) 

Thus Pi + \t) for any time t between and +oo can 
be approximately expressed by adiabatically connecting 
Eq. EH with (|BT7|l as follows, 



- +iA + Ae<*(t) 



A" 



(B18) 



Both Eqs. and (HIHI lead to the correct P a (oo) 

for steady states, 



P"(oo) = -- de 

7T . 



1 



e - h D (oo) +iA + Ae Q (oo) 



A a ,(B19) 



If the external applied voltage assumes a step-like form, 
for instance, AV a (t) = -Ae a (t) = AV a (l - e~^ a ) with 
a — > + , and hn(t) is not affected by the fluctuation 
of cT£)(i), Eq. (|B18|I would recover exactly Eq. (|B15|) . 
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In other cases, Eq. i|B18(l provides an accurate and ef- 
ficient approximation for Eq. (|B15|) . so long as AV a (t) 
do not vary dramatically in time. Since the integration 
over energy in Eq. (|B18|) can be performed readily by 
transforming the integrand into diagonal representation, 
Eq. (|B18p is evaluated much faster than Eq. (|B15|) . Due 
to its efficiency and accuracy, Eq. I|B18(1 is combined with 
Eqs. (jB9(l - HB12l) to form the WBL approximation for the 
dissipation functional Q^ BL , and thus recovers Eq. i|21|) 
ofSec.|V| 

As discussed in Sec.|V] Q^f BL (t) depends explicitly on 



AV a (t), cr D (t), h D (t) and A Q , where h D {t) is directly re- 
lated to Pd( v , t) by the Poisson equation on D subjected 
to boundary conditions AV a (t), and A Q are associated 
with the DOS of D near the surfaces S a . Therefore in 
practice Q^f BL is a functional of pd{y,i) only, i.e., 



Qa (t) = Qa <TD\pD],h D [a D \p D ],t], 



A a [p D ,t},AV a [ PD ,t},t 



(B20) 
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